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Abstract 

Constructing models to discover physics underlying magnanimous data is a traditional strategy in data 
mining which has been proved to be powerful and successful. In this work, a multi-optimized recurrent 
neural network (MRNN) is utilized to predict the dynamics of photosynthetic excitation energy transfer 
(EET) in a light-harvesting complex. The original data set produced by the master equation were trained 
to forecast the EET evolution. An agreement between our prediction and the theoretical deduction with an 
accuracy of over 99.26% is found, showing the validity of the proposed MRNN. A time-segment polynomial 
fitting multiplied by a unit step function results in a striking consistence with analytical formulations for 
the photosynthetic EET. The work sets up a precedent for accurate EET prediction from large data set 
by establishing analytical descriptions for physics hidden behind, through minimizing the processing cost 


during the evolution of week-coupling EET. 
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I. INTRODUCTION 


Many physical laws were discovered based on 
data, such as Kepler’s three laws via large da- 
ta of celestial body motion observed by Tycho. 
Understanding the temporal evolution of photo- 
synthetic excitation energy transfer (EET) in a 
light-harvesting complex is an important topic of 
broad interest due to its nearly 100% photosyn- 
thetic conversion efficiency, providing an ideal op- 
tion for mitigating energy crisis|1—3]. Exact nu- 
merical simulations of the dynamics of EET in 
a light-harvesting complex, on the other hand, 
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requires enormous computational resources|4, 5], 
which tends to grow exponentially with the num- 
ber of simulated time steps and system size. 
Even though many techniques, such as the hierar- 
chy of equations of motion (HEOM) technique|6, 
7], multi-configurational time-dependent Hartree 
(MCTDH)[8], stochastic Liouville-von Neuman- 
n equation|9], quasi-adiabatic propagator path- 
integral (QUAPI)|10], and path-integral Monte 
Carlo[11] are available, they are inappropriate for 
examining long-time quantum dynamical evolu- 
tion. 


The current state of a quantum system is most- 
ly determined by and rooted in its early stages 
of evolution, which enables us to learn long-time 


evolution of EET in a light-harvesting complex 
from short-time dynamics without the costly and 
direct long-time simulations. Once a memory k- 
ernel is acquired, the Nakajima-Zwanzig general- 
ized quantum master equation (GQME)[12] pro- 
vides a broad and formally reliable prescription 
for achieving this goal[{13]._ Nonetheless, solving 
the GQME and directly computing the memory 
kernel for an arbitrary system is still challeng- 
ing. The transfer tensor method (TTM) can solve 
the GQME somehow, but it requires an external 
numerical methodology to provide a set of dy- 
namical mappings|14—16]. The interplay between 
machine learning and quantum physics, on the 
other hand, altered the current situation[17—19] 
by providing new concepts for modeling the evo- 
lution of EET in a light-harvesting complex[20], 
i.e. intuitively and directly learning from data set 
without explicit theoretical construction. Artifi- 
cial neural networks (ANNs) [21] have shown that 
complex functional dependencies in time series 
can be learned directly from data[22, 23]. This 
evades the great efforts to make theoretical analy- 
ses, which sometimes become unjustified (e.g., the 
weak coupling limit) or difficult to derive if just 
based on phenomenological observations. ANNs 
have been demonstrated to be capable of solving 
the master equations governing the dynamics of 
long-time dissipative open quantum systems|24]. 


The recurrent neural network (RNN) in par- 
ticular has an exceptional capacity to interpret 
an intricate temporal behavior. It can preserve 
historical information for future purpose of pre- 
diction by building a feedback loop that receives 
both the current stage’s input and the previous 
step’s output|25]. However, concern of gradient 
vanishing or exploding behavior due to multiple 
iterating operations, on the other hand, limits the 
use of RNNs in long-time scale applications|26]. 
To address this flaw, this work employs multi- 
optimized recurrent neural networks (MRNNS) 
rather than long short-term memory recurren- 


t neural networks (LSTM-NN)[27, 28], to mod- 
el the long-term dependencies of the time-series 
data set, by storing the key information, and pre- 
dicting the future data that is not currently avail- 
able. MRNNs are also used as propagators of the 
time-dependent master equations to regulate the 
light-harvesting complex over a range of time s- 
cales. 

The rest of the paper is organized as follows. 
The quantum processes and the common mas- 
ter equations for photosystem II reaction cen- 
ter (PSII-RC) are described in Sec.2.A, and a 
sample RNN architecture with optimized hyper- 
parameters is introduced in Sec.2.B. Results are 
discussed in Sec.3, where we validate the learn- 
ing model in the interval [0. 80]fs (Sec.3.A) and 
predict the EET evolution process in the range 
[80, 500]fs (Sec.3.B), using the polynomial fitting 
and the analytical expression(Sec.3.C). Finally, 
conclusions and outlook of the this work are sum- 
marized in Sec.4. 


Il. THEORETICAL MODEL AND MULTI- 
OPTIMIZED RECURRENT NEURAL NET- 
WORKS 


A. Theoretical model for Photosystem II 
reaction center (PSII-RC) 


A typical photosystem II reaction center 
(PSH-RC) often seen in purple bacteria and 
oxygen-evolving organisms (cyanobacteria, al- 
gae, and higher plants) comprises six pigment 
molecules located at the core of the complex with 
two symmetric branches of protein matrix[29]. 
Four chlorophylls of them (special pair PD, 
PD, and accessory ChID,, ChlD2) and two pheo- 
phytins (PheD,, PheDg2), are parallel distributed 
in these two branches of protein matrix|30]. The 
pair of chlorophylls, PD;and PD», located at the 
center of the PSII RC act as the primary electron 
donors, forming two excited states denoted as 


FIG. 1. (Color online) Energy-level framework for the photosyn- 
thetic RC with two load-transitions |a;)—>|8;);-1,2. The electronic 
transitions from ground state |g) to two coupled dipoles |e1) and |e2) 
are induced by the high temperature photon bath, while the low tem- 
perature phonon bath drives charge transfer from |ez) to |a;)i=1,2, 


and |3;);=1,2 to |g) with a termination of electronic circulation. 


le12) in Fig.1. Two pheophytin pigments, PheD, 
and PheD, couple to the rest of the molecules 
and act as the electron acceptors|a,) and |ag), 
respectively|31], as shown in Fig.l. This is an 
energy-level framework abstracted from the PSII- 
RC[2], where |,)(|G2)) is a positively charged 
state after an electron is released and |g) is a 
ground state. 


After absorbing a solar photon, an electron is 
excited from |g) to |e2) with a transition rate yp, 
where the excited electron may transit to |e,) s- 
tate at a rate ye. Then the excited electron is 
transferred to the acceptors by emitting a phonon 
via two pathways: |e1)(\e2))— Jai) or |e1)(|e2))—> 
ag) at emission rates yic and Yz, respectively, 
with |a) and |az) being the ion-pair states in 
these two pathways. Such a process is accompa- 
nied by a spatial separation of positive and nega- 
tive charges induced by the release of the excited 


electrons to the plastoquinone molecule, leaving a 
hole in the dimer. Then the electrons are released 
from the two acceptors PheD, and PheD» denot- 
ed by |a1) > |61) (Path 1) and Jaz) > |82) (Path 
2) with respective rate |[;)@=1,2), providing ener- 
gies for possible work. Finally, the electron re- 
turns to the primary electron donor via |/5;)(|G2)) 
— |g). As an alternative pathway to the two- 
step processes |a)(|a2)) —|61)(|62)) —> |g), the 
acceptor-to-donor charge recombination can also 
be made by directly bringing the system back to 
the ground state |g) without producing curren- 
t with rates xT j212, where x is a dimensionless 
fraction|[2] describing the radiative recombination 
rate of the two pathways. 

It is significant to gather a data set, against 
which the learning model can validate to make 
prediction. In this work, the data used for train- 
ing is collected from the initial stage of above 
PSH-RC, which is generated by the weak contact 
of the system with the environment. 

The unitary evolution of the electron transfer 
can be described by an electronic Hamiltonian 
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via the Lindblad-type master equation 
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where the strength of the interaction with the en- 
vironment is comparable with the internal inter- 
actions inside the system. The Lindblad-type su- 
peroperators in Eq. (2) are listed as 
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Equation (3) describes the effect of the high 
temperature reservoirs with n, denoting its av- 
erage photon numbers. The low tempera- 
ture reservoir has the average phonon number- 
S Nic = [exp( #2) — 1|} in Eq. (4), where 
Yiic™ViclYjje=Yjc) are the spontaneous decay rates 
from level |e2) to level |a;}(i = 1,2) respective- 
ly. 
scribes the effect of Fano interference. 


The cross-coupling Yije with Yije=YVjie de- 
Similar- 
ly, another low temperature reservoir is described 
by Eq. (5) with na=lexp( 2”) — 1]! be 
ing the cold reservoir phonon numbers. Here 
Pije=T jic is defined as Pije = Mm y/T iT je describ- 
ing the Fano interference induced by the spon- 
taneous decay rates, Pge=Ise(T je=! je) G, j=1,2), 
from level |8; 1,2) to level |g) with ņı denoting 
the quantum interference robustness. In Eq. (6), 
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and Eq. (8). Next, 1 million training data set- 
s within 100 fs will be collected from the densi- 
ty matrix element equations (see the Appendix) 
with parameters|2] listed in Tab. (1). 


B. Multi-optimized recurrent neural net- 
works(MRNN) 


In this work, the dynamic evolution of popu- 
lations on each energy level in Fig. 1 was utilized 
to illustrate EET evolution, and the long-term 
evolution of EET in the proposed PSII-RC will 
be predicted by a recurrent neural network (RN- 
N) tuned by four hyperparameters. The simple 
RNN[21, 32] is a type of neural network that is de- 
signed to learn data sequences such as time series, 
as illustrated in Fig. 2. To further understand 
how they function, it is reasonable to compare 
them to more traditional feedforward neural net- 
works. In feedforward neural networks, the input 


TABLE I. Model parameters used in the numerical calculations. 


Values Units 
Ee, — Eg = Ee, — Eg 0.185 eV 
Ee — Ea, = Ee. — Eog 0.2 eV 
Eg, = Ey = Eg, = Eg 0.2 eV 
Yh 2.48 x 107 eV 
Ye 0.025 eV 
Yle = Y2e 0.012 eV 
rı =I, 0.124 eV 
Tye = Tae 0.0248 eV 
Ta 0.026 eV 
Nh 6000 
Ne 0.46 
x 0.2 
m 0.25 


data set is propagated step by step via numer- 
ous intermediary layers, and training is accom- 
plished by updating the weight matrices as well as 
the vectors before reaching the final output layer, 
such that the neural network learns certain de- 
sirable input-output correlations. In feedforward 
neural networks, the input data set is propagated 
step by step through multiple intermediary levels, 
and training is performed by updating the weight 
matrices and vectors before reaching the final out- 
put layer, so that the neural network learns some 
desired input-output relationships behind the da- 
ta. A feedforward network might theoretically be 
used to handle temporal data. 


However, a feedforward neural network is 
hardly the best solution since the number of free 
parameters rapidly grows with the number of time 
steps. As an alternative option, RNNs handle this 
issue with a cyclic connection design in which the 
update rule for the hidden layers at time t is de- 
cided not only by the current states S;, but also 
by the states at previous times(t — 1). t— 1, t, 
t+ 1 are the time series, £41, Zt, 441 represent 


FIG. 2. 


work(RNN) model. 


The architecture of the simple recurrent neural net- 


the input sample data, and S;_1, St, St+1 repre- 
sent the memory of the input sample x; at time 
ti, and O1, O; and O1, respectively, represent 
the status information stored at the current time. 
The RNN update rule is given by 


Si = fUr + W Si), 
O: = g(V S). 


(9) 
(10) 


When the data set propagates forward in Fig. 2, 
the three weight matrices W (input weight), U 
(input sample weight), and V (output sample 
weight) do not depend on time t, i.e., weight 
sharing. The activation functions f and g then 
take into consideration the free parameters W, U 
and V, as well as previous memory. The acti- 
vation function f can be tanh, relu, sigmoid or 
other functions, while g is usually softmax and 
other functions. Therefore, the state S;_; of the 
previous period will participate in the prediction 
of the state S;. In other words, the input data 
and the hidden state of the previous time step 
will be calculated using the weight matrices so as 
to generate the hidden state of the current time 
step[33, 34]. 

Because of these natural preferences, the 
Markovian assumption of independent data 
points at multiple time steps is problematic in 
the application of RNN. However, if the machine 
learning model favors noise and minutiae, and ig- 
nores the general trends and patterns in training 


the data set, it will perform well on trained da- 
ta but poorly on new data, a phenomenon termed 
overfitting. In this work, we will use a simple RN- 
N with a multi-optimizer to predict the quantum 
evolutions of EET in light-harvesting complexes. 
Additionally, the multi-optimized hyperparame- 
ters are addressed further below. 


1. Learning rate regulator (LRR) 


A learning rate regulator (LRR) will be fit- 
ted to the simple RNN, which will determine how 
much the model parameters are updated in each 
iteration and how far down the gradient the pa- 
rameter moves with each update[35, 36]. The ex- 
ponential attenuation regulator, cosine LRR, pre- 
heating regulator, and LR attenuation regulator 
are all commonly used LRRs. A successful LRR 
should be able to optimize the model while avoid- 
ing overfitting and underfitting. If the LRR is set 
too high or too low, the model may diverge or 
converge slowly, requiring more training rounds 
to get the optimal solution. 


To simplify the task, we frequently define an 
LRR at the starting stage of training. During 
training, the model computes the gradient of the 
loss function to determine the updating direction 
of the parameters, whereas the LRR dynamically 
updates LRR. The performance of the model is e- 
valuated using both the training and verification 
sets. Because the time series for the quantum evo- 
lutions of EET are projected to be long, a faster 
decay rate LRR, i.e., the exponential decay LRR, 


—epoch 
10 


LRR = 0.001 x exp( ), (11) 


will be employed to this RNN. Here the epoch 
parameter means the iterations during the train- 
ing. 


2. Early stop function 


In addition to the previously stated anti- 
overfitting strategies, an early stop function|37] is 
employed to decrease the overfitting. As implied 
by its name, the early stop technique completes 
training before the algorithm overfits and obtains 
the optimal global outcome, resulting in robust 


generalization performance, as shown by 


Eopt(t) := Miny <= Eyal), (12) 
GL(t) = 100 x poe —1], (13) 


where Eop(t) is the ideal verification error set 
as a function of the number of repetitions t, and 
GL(t) is the generalization loss evaluating the 
rate at which the generalization error grows in 
comparison to the previous lowest error. When 
the generalization error is large, an early stop is 
preferable since it indicates that the model has 
been fitted. Such an ending of training is judged 
by a threshold of GL(t). 
nique’s halting criterion is classified into three 


The early stop tech- 


types: the first, second, and third. In this work, 
the first type of stop rule is employed to determine 
the loss function and accuracy on the verification 
set. 


3. Regularization 


Regularization with the additional penalty 
terms into the model’s loss function is becoming 
a popular strategy for reducing model overfitting 
and improving model generalization in machine 
learning. There are two typical methods of reg- 
ularization, namely L1 and L2([38, 39], evaluated 
by two weight parameters w* and w following and 
before the update, respectively. They are given 


by 
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where the sum of the absolute values and the 


FIG. 3. Recurrent Neural Network without Dropout the structure 
diagram of standard neural network (a) and its modification with 


dropout (b). 


square of the parameters have been added to the 
loss function, respectively[39]. The first terms on 
the right side of the two equations represent the 
original loss functions, the second terms represent 
the regularization parameters A, and the summa- 
tion is carried out over the numerous model own- 
ership parameters. L1 regularization (14) causes 
the model to create sparse parameters (some pa- 
rameters are set to zero), whereas L2 regulariza- 
tion (15) causes the model parameters to be near 
to but not completely zero. As a result, L1 reg- 
ularization creates sparse solutions and is more 
suitable for feature selection, whereas L2 regular- 
ization produces a greater number of tiny non- 
zero weights. In this study, L2 regularization will 
be used to prevent overfitting in the long-term 
training data set. 


TABLE II. Hyperparameters used in the RNNs train- 


ing 


L2 LR epochs Dropout 
Para, 0.01 0.0001 139 0.11747474747474748 
Pasay 0.01 0.0001 140 0.09779519760654846 
Pee, 9.01 0.0000001 155 0.32499999999999996 
Pese. 9.01 0.0000001 125 0.32499999999999996 
Pop 9.01 0.0000001 122 0.35 


4. Dropout 

In this work, the Dropout regularization 
technique[40—42] is also applied to the simple RN- 
N. The input training data is carried forward 
through the neural network, the estimated loss is 
propagated back, and the parameters are changed 
using the gradient descent approach, as shown in 
Fig. 3(a). In Fig. 3(b), the Dropout parameter P 
is used to deactivate certain neurons with a spe- 
cific probability during the forward propagation, 
and the parameters are updated using the gradi- 
ent descent approach. This method is repeated 
multiple times so as to properly alleviate overfit- 
ting during training in this proposed MRNN. 


II. RESULTS AND DISCUSSIONS 
A. Training the multi-optimized recur- 
rent neural network (MRNN) 


By progressively incorporating a few hyperpa- 
rameters, such as early stop function, L2 regu- 
larization method, LRR, Dropout, and Bayesian 
optimizer, a multi-optimized RNN model is con- 
structed, which, as a distinguishing feature, uti- 
lize the Lindblad-type master equation Eq. (2) 
to gather data from the proposed PSII-RC, along 
with certain model parameters listed in Tab. (1). 
We used Eq. (2) to produce a data set of 1 million 


data points in 100 fs, which are divided into train- 
ing and test sets at 4:1 ratio. The first 800,000 
data points serve as the training set, while the 
remaining 200 000 act as the test set, using the 
hyperparameters listed in Tab. (II). Fig. 4(a) dis- 
plays the evolution of the training set fed in- 
to this proposed multi-optimized RNN learning 
model (Original codes in SM1) during the inter- 
val [0, 80]fs. 


FIG. 4. (Color online) (a) Evolutions of the training population 
at all excited states within the interval[0, 80]fs. (b) A comparison 
of learning quality between our proposed learning model and the 
collected testing set over the range [80,100]fs. The hyperparameters 
listed in Tab.(II) and the original codes in SM1, SM2 and SM3 are 


used for the calculation. 


To assess the validity of the proposed learn- 
ing model, we predict the evolution of EET (solid 
curves) over a time range from 80 to 100 fs with a 
comparison to the test set (dotted curves) collect- 
ed from the aforementioned PSII-RC within the 
same range, indicating a good agreement between 
them, as illustrated in Fig. 4(b). This ensures the 
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high accuracy of this learning model composed 
of the simple RNN and some optimizers(Original 
codes in SM1). Supported by these results, the 
proposed multi-optimized RNN learning model is 
expected to be able to anticipate the evolution of 
EET from 80 to 500 fs. 


B. EET Predicted by the MRNN 


In Fig. 5, the vertical dashed black lines at 
the 80th fs shows the temporal starting point of 
our prediction achieved by gradually incorporat- 
ing hyperparameters into calculation. The cut-off 
point of 80 fs can be used to evaluate the accuracy 
of the predictions by checking how precisely the 
data generated by the prediction coincide with 
the training data set. After examining Figs. 5(a)- 
(e), the evolutionary history of each population 
indicates that our predictions vary against hy- 
perparameters in the following periods. L2 reg- 
ularization, on the other hand, has little effect 
on EET prediction, as indicated by the nearly i- 
When all the 
hyperparameters are added into the simple RNN, 


dentical curves in Figs. 5(a)-(e). 


there is a notable convergence between the train- 
ing and prediction values at 80 fs, corresponding 
to the red curves in Figs. 5(a)-(e)(Original codes 
in SM2). 

Consider the prediction of EET on the excited 
state |e) in Fig. 5(a), where the hyperparame- 
ters and optimizers are added one by one with- 
in the RNN architecture. Three layers with 128, 
64, and 32 neurons, respectively, are the essential 
characteristic neural network architecture. In the 
absence of optimizers and hyperparameters, the 
simple RNN timing prediction model fails to yield 
significant physical results, as shown by the curve 
“ Simple RNN ” in the insets of Fig. 5(a), as ev- 
idenced by the almost coincided curve with the 
curve “Simple RNN ”. The time forecast provid- 
ed by the later added early stop function does 
not satisfy, which can be drawn by the almost 
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FIG. 5. 


shown in (a) to (e) via the multi-optimized RNN in the interval [80, 


(Color online) Prediction of the evolutive EET were 


500]fs, and the predictive and theoretical (within [0, 500] fs) EET 
evolutions were collected in (f) with identical parameters to those in 


Fig.4. 


coinciding curve with the curve “Simple RNN ”. 
Although the apparent physical rule is predict- 
ed, their prediction accuracy is practically iden- 
tical, as shown by the green and purple curves 
in the inner inset. The evolution curve of |e1) 
tends to level off when the Dropout parameter 
(P=0.32499999999999996 is applied for further 
optimization, as demonstrated by the red curve. 
The overlapping purple and green curves in the 
insets of Fig. 5(a) indicates that the development 
of EET is not sensitive to the L2 regularization 
on the excited state |e;). 

In Fig. 5(c), the red curve achieves per- 
fect docking with the 
80 fs with the Dropout parameter set as 
P=0.11747474747474748 and the initial LR= 
0.001. As it takes time to manually adjust the 
hyperparameters of the neural network when 


training data at 


optimizing the prediction of the dynamical 
population on jaz), the Bayesian optimizer is 
implemented into the neural network to find the 
ideal number of layers and neurons. Finally, a 
neural network composed of 6.23223034545932 
layers with 124.19253861421566 neurons per 
layer was employed, corresponding to the red 
curve shown in Fig. 5(d). The same neural 
network architecture parameters are employed 
to forecast |b), and the inset in Fig. 5(e) clearly 
shows the roles each hyperparameter plays in 
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FIG. 6. 


by the population sum in the interval [80, 500]fs. 


(Color online) Evolution of prediction accuracy defined 


To validate the predictive accuracy of the pro- 
posed MRNN in this work, the evolution behav- 
iors of EET are theoretically simulated by using 
Eq. (2) in the temporal range of [0, 500] fs, as 
shown by the dotted curves in Fig. 5(f), while 
the optimal prediction results are re-plotted as 
curves with points in the time range 80 to 500 fs, 
as shown by the red curves in Figs. 5(a)-(e) for 
comparison(Original codes in SM1). 

It is found that, as illustrated in Fig. 5(f), the 
theoretical calculation and the curves predicted 
by MRNN for photosynthetic EET perfectly co- 
incide, indicating the validity of our proposed M- 
RNN. In addition, the predicted values of each 
population were added together and divided by 
one during [80,500 fs], to evaluate EET predictive 
precision for the system due to their theoretical 
sum of one unit. Although the predictive pre- 
cision decreases against time for this photosyn- 
thetic model, it reaches a low of 0.9926 at 500fs, 
exhibited by an arrow in Fig. 6(Original codes in 
SM1). Nevertheless, this accuracy outperforms 
most RNN learning models, such as those report- 
ed by Refs. [43, 44]. At the same time, because 
our MRNN learning model is born from the sim- 
ple RNN, long-term prediction inevitably will re- 
veal the inherent flaw of short-term prediction, re- 
sulting in a progressive drop in accuracy, as seen 
by the decreasing feature of Fig. 6. 


C. 
pression for the predictive EET 


Polynomial fitting and analytical ex- 


The purpose of this work is to discover the 
physical rules underlying photosynthetic EET da- 
ta. Fitting the predicted findings based on da- 
ta with a polynomial is an effective mathemat- 
ical tool to get an analytical understanding of 


4 =» Predictive EET evolution Pee, 
1 = += Polynomial fitting pee: 


200 


ered _ 
eee ve 


i —. Predictive EET evolution pee, 
== Polynomial fitting Peze 


200 400 


l —. Predictive EET evolution Paa, 
== Polynomial fitting Pao, 


—+ Predictive EET evolution Paa, 
== Polynomial fitting Paza, 


4 —- Predictive EET evolution ppp 
\ === Polynomial fitting po» 


(Color online) Polynomial fitting (red dotted 


EET. Such a fitting usually relies on the least FIG. 7. 
squares method and lowering the error sum of  curves)comparised with the predictive evolution (blue dash-dotted 
squares[45, 46]. In Python, the Polynomial Fea- 
tures function is a utility in the scikit-learn library 


curves) of the evolutive EET within the interval of [80, 500]fs. 
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to generate the polynomial features so as to form 
a new feature matrix. Finally, the feature matrix 
is returned for further model training and fitting. 
When a polynomial of order 2 is provided for a 
one-dimensional feature A, Polynomial Features 
generate a new feature matrix containing A to 
the first and second powers. If the original fea- 
ture has many dimensions, the resulting feature 


matrix will have power combinations for each of 
them. 

Figure 7 displays the differences between the 
predicted curves and the fitting curves using Poly- 
nomial Features in scikit-learn library. The maxi- 
mum degrees of polynomial fitting for each energy 
state all are 5 from (a) to (e) in Fig. 7 (Original 
codes in SM3). And the analytical polynomial 
fittings to Figs. 7(a)-(e) are also given by 


Pee, (t) = —0.0022184223655202433¢ + 1.3936336816327624 x 10~*¢? 


Peres (t) = 


—4.3658376701560287 x 1078t? + 6.694506984131638 x 107141 
—4,.008815848721703 x 107145 + 0.33721254828307511, 
0.0055891350945578t — 3.50638450439065 x 1075: 


+9.953236151938327 x 107-4? + 0.24964385363485108, 


Poara, (t) =4.317228690621229 x 107° — 3.0372753727222955 x 107174? 


+1.0249801765297704 x 10~7°t° + 0.00000002169596502, 


Pose (t) = —1.8447722039325343 x 107% + 1.0562606143359824 x 10782? 


—2.4712186236144212 x 10~1’¢° + 0.00021705673079966, 
—0.003319998227995093t + 2.0476938571365984 x 107 °t? 


Poo (t) = 


(16) 
+1.091124267760133 x 10~t? — 1.6659867620310589 x 1071044 

(17) 
+1.0210132754732106 x 107'#3 — 1.647224363731816 x 1071"t4 

(18) 
—3.040707800147938 x 1071443 + 4.356557804662399 x 107114 

(19) 
—6.313410066086038 x 1078+? + 9.606627816067419 x 1071144 

(20) 


—5.735689700969715 x 107 ‘44° + 0.41192131944573884. 


The polynomial fitting curves of Pe,e,, Peres Pasar: 
and pæ agree with the predicted curves fairly well, 
as shown in Fig. 7. In contrast, the fitting re- 
sult of Paia, is rather poor, as shown in Fig. 7(c). 
Overfitting is visible in its dynamic population 
progress. Therefore, an alternative fitting tech- 
nique should be found for paio- 

Because the characteristic features of an RN- 
N lies in its internal (hidden) loop memory, it 
is understandable that a dynamic state contain- 
s information on all previous input as it evolves 
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through the data sequence, as shown in Fig. 2. As 
the polynomial order is increased, the polynomial 
features over adapt to the past information but 
perform poorly on fresh data. This enlarges the 
changing amplitude and change rate that lead to 
overfitting. An effective method is to derive in- 
formation based on all previous inputs. For that 
purpose, a time-segment strategy is implemented, 
i.e., polynomial fittings in different time ranges 


—- Predictive EET evolution pa,a, 
+ Polynomial fitting Paa, 


Pe a a 
1 


FIG. 8. (Color online) A comparison between the fittings by the 
time-segment polynomial multiplied by unit step function (red dot- 
ted curves) and the predictive evolutions (blue dash-dotted curves) 
in [80, 500]fs. Time-segment polynomial fittings in comparison with 
the predictive evolutions shown by the insets(a1) in [80, 250]fs and 


(a2) in [250, 500]fs. 


are carried out and a unit step function 


t < 250 0 
t) = 
t >= 250 Alt) f 


is incorporated to the final polynomial expression 
in each time-segment. The fitting time is divid- 
ed into two intervals [80, 250] fs and [250, 500] fs, 
and the polynomial fitting in each interval is mul- 
tiplied by a unit step function given by Eq. (21). 
Their final fitting is the sum of the two polyno- 
mial fittings, and given by 


Paia (t80-500) = Paia: (ts0-250) * f(t) 
+para: (t250-500) * fa(t). (22) 


The curves in Fig. 8 shows the fitting result- 
s for the population on the state |a) using the 
aforementioned procedure (Original codes in S- 
M3). The insets Figs. 8(al) and (a2), show a 
comparison between polynomial fitting and pre- 
dictions in the intervals [80,250]fs and [250,500] 
fs. The highest degrees of the polynomial fitting 
is set as 5 in Figs. 8(al) and (a2). The over- 
lapping between the dotted lines (red) and dash- 
dotted lines (blue) demonstrates the precision of 
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this polynomial fitting at various intervals. Fur- 
thermore, the analytical functions derived by the 
time-segment polynomial fitting are, respectively, 
given by 


Pare, (tgo—250) = 1.7388382452293818 x 107°t 
—1.8538908793974676 * 107"? 
+9.883817236997048 x 1071" 
—2.6256355830257586 x 107161“ 
+2.7737661365102853 « 1071°{23) 

Para, (t250-500)= 2.500773638209706 x 107 "t 
—1.2744394675654058 « 10742? 
+3.227678944102217 « 10717’ 
—4,0616638358969016 * 10701“ 
+2.031479651834943 * 10~”t124) 


The comparison made between the whole pe- 
riod fitting and predictive evolutive curves in 


< “P2180. 500|fs for the dynamics populations on the 
> — 


states |a,) demonstrates that the time-segment 
polynomial fitting can well overcome the overfit- 
ting in Fig. 7(c). 

Figure (9) depicts the total fitting loss rate 
compared to the predictive results, a physical 
quantity assessing the precision of the fitting tech- 
nique employed in this work. The loss rate curve 
in Fig. (9) shows an oscillating behavior in the 
time interval of interest. Although a jumpy loss 
rate can be seen both at the initial and the final 
stages of the time interval, it still remains around 
an order of 1075, ensuring a high accuracy for this 
MRNN prediction utilized in this work. 


IV. CONCLUSION AND OUTLOOK 


In summary, fed by the original PSH-RC data 
set, a MRNN strategy is proposed to forecast the 
EET evolution, with an accuracy of over 99.26% 
within 500fs if compared to the theoretical deduc- 
tion. The polynomial fitting is also implemented 
for the EET evolutions so as to get analytical re- 
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FIG. 9. (Color online) Total loss rate of the polynomial fitting 


versus the prediction within the range of [80, 500]fs. 


sults. The predicted EET evolutions are also sub- 
jected to time-segment polynomial fitting multi- 
plied by a unit step function, and the analyti- 
cal formulations with high precision demonstrate 
the closeness to physical law. The results reveal 
that the proposed MRNN is a valid and power- 
ful data mining tool for forecasting the evolution 
of EET in a light-harvesting complex, and this 
study establishes a precedent for predicting EET 
from data using MRNN. A comparison to exper- 
imental data in the future is expected to assess 
| 


Pee = Yelne + 1)pere 


the validity of this learning model. 
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VI. APPENDIX 


The density matrix dynamic element equations 
are given by 


— aces = Plina Eg DD) sees ~~ MnPgg], 


Peres = Yelne + L)Peze: = Nepezez] E Viel(N1¢ =F 1) Bees = Nei Paran] 
—Y2e|(Me + 1)Pezez ca Nie bands -+ 2yi2cNce1 RelPoiazl, 


Porai = Viel Me + 1)Pezes = Mines | ~ Fiche Pas | (1+ A) eis 

Pozo. = Viel(M1e + 1)Pezez — NicPoror] — VrzeM1eRe[Para2] — (1 + AYP 2P ares, 

baias = 1A tPaies 9 (Ne + V2c)NicParaz + i2e[2(Mre + 1)pesez ~ NicPazaz — MicPoara] 

Ppp, = FiPoaro — Vrel(M2e + 1)P6,8, — N2cPgg] — Fizn: + 1) Reloge], 

P86 = T2Paza2 — V2e|(M2e + 1)P8263 — NzcPpog] — Vi2e(M2e + 1)Relps, 65], 

Ppp. = IAP — “Oe + Pze) (nze + 1) P6162 — SP rl(r + 1)ppia, + (Nze +1) p26. — 2N2cP99) 
Pog = L= Perei T Pezes — Paroi T Paras T Phib T PB2B21 


where Ay = Ea — Ea, and Aa = Eg, — Eg, are the 
splitting of the states |a,)(\a2)) and |81) (|820). 
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We utilize the equations to simulate dynamics of 


EET in PSII-RC. 
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